Volume 112, Number 2, March-April 2007 

Journal of Research of the National Institute of Standards and Technology 

[J. Res. Natl. Inst. Stand. Technol. 112, 75-94 (2007)] 



He Thermophysical Properties: New Ab Initio 

Calculations 



Volume 112 



Number 2 



March-April 2007 



John J. Hurly 

National Institute of Standards 
and Technology, 
Gaithersburg, MD 20899-8360 

and 

James B. Mehl 

K T Consulting, Inc., 

P.O. Box 307, Orcas, WA 98280 

john.liurly@nist.gov 
jmehl@rockisland.com 



Since 2000, atomic physicists have 
reduced the uncertainty of the helium-heh- 
um "ab initio" potential; for example, from 
approximately 0.6 % to 0.1 % at 4 bohr, 
and from 0.8 % to 0.1 % at 5.6 bohr. These 
results led us to: (1) construct a new inter- 
atomic potential 0q7, (2) recalculate values 
of the second virial coefficient, the viscos- 
ity, and the thermal conductivity of He 
from 1 K to 10,000 K, and (3), analyze the 
uncertainties of the thermophysical proper- 
ties that propagate from the uncertainty of 
0Q7 and from the Bom-Oppenheimer 
approximation of the electron-nucleon 
quantum mechanical system. We correct 
minor errors in a previous publication [J. 
J. Hurly and M. R. Moldover, J. Res. Nat. 



Inst. Standards Technol. 105, 667 (2000)] 
and compare our results with selected data 
pubUshed after 2000. The ab initio results 
tabulated here can serve as standards for 
the measurement of thermophysical prop- 
erties. 
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1. Introduction 

In 2000, Hurly and Moldover published a compre- 
hensive report on the application of fundamental 
physics to the calculation of the thermophysical proper- 
ties of low-density helium [1]. The present paper is an 
extension to and update of parts of that paper. We 
developed a new model potential for the interaction of 
helium atoms, (f>Q-]{r\ based on the most recent theoret- 
ical values of <^r). This potential was used to calculate 
several important properties of "^He: The density virial 
coefficient B{T) and its first two temperature deriva- 
tives, the zero-density viscosity, and the zero-density 
thermal conductivity 

Our improved potential and calculations have signif- 
icantly reduced the uncertainty of the thermophysical 
properties of helium. For example, at 300 K, the uncer- 
tainty of the second virial coefficient is now 1/7 of that 



reported in Ref [1] and the uncertainty of the thermal 
conductivity is 1/3 of that reported in Ref [1]. 

The new potential includes the diagonal correction to 
the Bom-Oppenheimer model (DBOC). In addition to 
the use of this correction, recent discussions of the adi- 
abatic model [2,3] recommend the use of atomic, rather 
than nuclear, masses in the calculations of atomic inter- 
actions. We have examined the sensitivity of the ther- 
mophysical properties to this replacement, as well as to 
the DBOC and to the uncertainties of the theoretical 
calculations of (p. 

In the temperature range 3 K to 933 K, helium gas 
thermometry [4] played a leading role in the formation 
of the internationally accepted temperature scale, ITS- 
90. Subsequently, improved gas thermometry has 
measured T- T^q, the differences between the thermo- 
dynamic temperature and ITS-90. Thus improved gas 
thermometry [5] may lead to a future, improved tem- 
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perature scale. Each form of gas thermometry (constant 
volume, dielectric, acoustic) requires the extrapolation 
of measured gas properties to zero pressure, where 
gases become "ideal." In this work, we use fundamen- 
tal principles to calculate the second density virial coef- 
ficient of helium B{T) and the second acoustic virial 
coefficient of helium ^aiX) with smaller uncertainties 
than can be achieved by direct measurements. Our tab- 
ulated values for B{T) and ^8^(7) can be used to con- 
strain the extrapolations to zero pressure; thereby lead- 
ing to more accurate values of the thermodynamic tem- 
perature. Acoustic gas thermometry also requires accu- 
rate values of the thermal conductivity, which we have 
tabulated for helium. Recently, May et al [6] have 
shown how to combine ab initio values of the viscosity 
of helium with comparatively simple viscosity-ratio 
measurements to obtain values of the thermal conduc- 
tivity of argon that are more accurate than can be 
achieved by direct measurements. Thus, our tabulation 
of the viscosity of helium will also facilitate more accu- 
rate argon-based acoustic thermometry. Finally, we 
mention programs to redetermine the Boltzmann con- 
stant [7] and to develop an atomic standard of pressure 
[8] based on accurate measurements of the dielectric 
constant of helium at the temperature of the triple point 
of water (7Vpw = 273.16 K). Both of these programs 
will benefit from the reduced uncertainties of 5(7). 

In the following sections, we first describe the poten- 
tial and the way it was developed. We summarize the 
quantum-statistical formulas used for calculating the 
thermophysical properties of interest, then describe the 
numerical procedures used for the calculations. We 
conclude with some comparisons of our theoretical 
thermophysical properties with recent experimental 
results. 

Standard notation conventions are followed in this 
paper. All quantum-mechanical formalism is expressed 
in atomic units except when noted otherwise. 
Interaction potentials are expressed in hartrees in the 
formalism, but converted to temperature units (K) for 
comparison with relevant literature. The CODATA- 



2002 values of the fundamental constants [9] were used 
in all calculations. 



2, Model Potential ^o? 

The potential model is expressed as the sum of a 
repulsive term and an attractive term 



<^07W^ 



kepW + ^attW, r^ <r<oo 



<^iepW = ^exp Xa,r" 



0attW = -^S 



/2«WQ, 



{5rf_ 

[ k\ 



-dr 



(1) 



(2) 



• (3) 



In these equations the cut-off radius Tq = 0.3 bohr is 
chosen to exclude the unphysical behavior of the poten- 
tial model at small r\ A= \ hartree {E^ defines the 
units; the a„ and 5 are fit parameters; the C2„ are fixed 
parameters; and the functions ^„ account for relativistic 
retardation. The attractive part of the potential is the 
sum of multipole attractive terms multiplied by the uni- 
versal damping functions of Tang and Toennies [10]. 

The dipole-dipole and higher multipole parameters 
Q for 77 < 10 (Table 1) are fixed at the values calculat- 
ed by Zhang et al [11]. The coefficients Q for helium 
with the mass of "^He were used in all calculations 
except those which investigated corrections to the 
Bom-Oppenheimer model. The coefficients C2„ for n > 
5 were estimated using the three-term recursion formu- 
la of Thakkar [12]. 

Zhang et al [11] include an extensive tabulation of 
previous calculations for comparison. The fractional 
differences between the fixed-nucleon parameters of 
Zhang et al and those of Bishop and Pipin [13] are 



Table 1. Attractive interaction coefficients [11] for helium atoms with ^e and infinite mass 
nucleii 



^e 



"He 



Q (hartree-bohr^ 
Cg (hartree-bohr^) 
Cio (hartree-bohr^^) 
Ci2 (hartree-bohr^^) 
Ci4 (hartree-bohr^"^) 
Ci6 (hartree-bohr^^) 



1.462122853192 

14.12578806 

183.781468 

3267.13274 

76501.2887 

2277412.86 



1.460977837725 

14.11785737 

183.691075 

3265.256092 

76571.26764 

2276292.717 
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2.6 X 10"^ for Q, 1.7 x 10"^ for Cg, and -9.5 x 10"^ for 
CiQ. If these fractional differences are taken as esti- 
mates of uncertainties, the total uncertainty in the 
potential is less than 3 x 10"^ K, and the total fractional 
uncertainty is less than 5x1 0"^, for r > 1 bohr. 

A fiirther, and more significant, source of uncertain- 
ty is the extrapolation formula used to estimate C2„ for 
n> 5 from the lower-n values of C2„. Thakkar [12] rec- 
ommends the use of either his Eqs. (29) or (33), with 
the latter more appropriate for helium (based on the 
value of QCio/Cg ). With the alternative formula, the 
estimated values of C12, C14, and Q^ are 1.3 %, 5.2 %, 
and 1 3 % larger. If these differences are used as esti- 
mates of the uncertainties of the corresponding poten- 
tial contributions, the total uncertainty in the potential 
is less than 4 x 10"^ K, and the total fractional uncer- 
tainty is less than 6x1 0~^, for r > 10 bohr. 

In principle, the use of the Tang-Toennies damping 
terms [10] is an additional source of uncertainty. 
However, these functions differ from unity only for r ~ 
10 bohr and below, where the quality of the fit potential 
can be judged directly by comparison with theoretical 
potentials. (See Fig. 2.) 

The retardation functions y^^/g, and/10 have been cal- 
culated by Chen and Chung [14]. Their results for/ are 
in excellent agreement with the calculations of 
Jamieson et al [15], whose results differ from those of 
Chen and Chung by a maximum fraction 1.5 x 10"^. 
The retardation fimctions satisfy /^(O) =1;/ decreases 
to V2 for r ^ 500 and approaches 328.47/r for large r;/ 
decreases to V2 for r ^ 660 and approaches 420.62/r for 
large r; and/o decreases to Vi for r ^ 8 10 and approach- 
es 508.43/r for large r. The functions /„ have the effect, 
for example, of converting the dipole-dipole interaction 
from a 1// dependence to a 1/r^ dependence. 
Retardation has, at most, a marginal effect on all terms 
except the dipole-dipole term. At r = 660 bohr, the ratio 
Cg/g/r^ to 007 is less than 3 x 10"^; similarly, at r= 810 
bohr, the ratio of Qo/o/r^^ to 0o7 is less than 4 x 10"^^ 
Accordingly, the factors /12, /u, and /g, were safely 
approximated as unity. The code for computing the 
potential uses cubic spline interpolation of the results of 
Chen and Chung [14] for/,/, and/o. 



The parameters 5 and Up -2 <j < 2 were determined 
by fitting the potential model (l)-(3) to selected theo- 
retical values weighted to account for their estimated 
uncertainties. The retardation functions /„ were set to 
unity in these fits. Several fits were made. The first was 
to the selected data set described below, and will be 
referred to as 0o7- The second and third fits accounted 
for the uncertainties of the theoretical values, also 
described below. The corresponding potentials are des- 
ignated 0o7±- An additional fit was made to the potential 
values without applying the diagonal Born- 
Oppenheimer correction [16]. The corresponding 
potential is 0nboc- The values of 5 and the Uj determined 
by these fits are listed in Table 2. 

2.1 Theoretical Values of ^ 

Table 3 lists values of the potential (p and their uncer- 
tainties based on our review of the recent literature 
[17-26]. The values selected for determination of 0o7(^) 
represent a compromise based on availability of calcu- 
lations at each r, the uncertainty claimed by the authors, 
and the internal agreement of various calculations for 
nearby r. The theoretical values were obtained within 
the Bom-Oppenheimer model for fixed nuclear separa- 
tions. Uncertainties were assigned to each of the select- 
ed values. When only a single datum was available, the 
authors' uncertainty estimate was used, provided that it 
was consistent with neighboring values; otherwise the 
uncertainty was adjusted upward. When several values 
were available at an r-value, generally the unweighted 
mean and standard deviation of the more recent calcu- 
lations was used. The upper-bound potentials of 
Komasa [19] were used only at small r, where they are 
in excellent agreement with the quantum-Monte-Carlo 
calculations of Ceperley and Partridge [17], which have 
much larger uncertainties. 

The diagonal Bom-Oppenheimer correction calcula- 
tions of Komasa, Cencek and Rychlewski [ 1 6] were 
interpolated using a cubic spline and added to the fixed- 
nucleon potentials. 

Relativistic [27] (+15.4 mK) and radiative [28] (-1.3 
mK) corrections to the potential have recently been 



Table 2. Variable (fit) potential coefiicients 



Potential 



(bohr^) 



(bohr) 



(-) 



(bohr"') 



^2 
(bohr"T 



3 
(bohr"^) 



007 


0.081212 


-0.28755 


2.14735 


-1.97272 


-0.051787 


1.992657 


007- 


0.097486 


-0.32441 


2.17654 


-1.98206 


-0.050505 


2.006175 


007+ 


0.065002 


-0.25089 


2.11837 


-1.96343 


-0.053050 


1.980020 


0nboc 


0.072490 


-0.26814 


2.13133 


-1.96754 


-0.052524 


1.985551 
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Table 3. Theoretical potential values used to determine 0o7- The model potential 0o7 was fit to the sum 
of the theoretical potential (j) and the diagonal Bom-Oppenheimer correction A(^boc ^i^ weighting 
equal to the inverse square of the uncertainty t/(0). When a single source is listed, the uncertainty is gen- 
erally that stated in the source. When multiple sources are cited, the unweighted mean and standard devi- 
ation of the set is used. In some cases, indicated by an asterisk, the uncertainty was adjusted upward to 
account for disagreement with neighboring values. 



r (bohr) 


0(K) 


t/(0) (K) 


^<fcBOC (K^) 


Source(s) 


1 


286435 


25 


158 


[19] 


1.5 


104320 


20 


36 


[19] 


2 


36144.6 


10 


11.8 


[19] 


2.5 


11962.0 


1.0 


4.1 


[19,23]* 


3 


3786.0 


20 


1.37 


[17,19,23,24] 


3.5 


1111.0 


1.0 


0.41 


[19,20,23]* 


4 


292.64 


0.10 


0.10 


[20-23,25,26] 


4.5 


58.400 


0.10 


0.009 


[19,20,23,24] 


5 


-0.500 


0.10 


-0.013 


[19,20,23,24] 


5.1 


^.534 


0.025 


-0.014 


[20,23] 


5.6 


-10.991 


0.011 


-0.012 


[20-26] 


6 


-9.671 


0.009 


-0.011 


[20,23]* 


6.5 


-6.887 


0.005 


-0.008 


[23] 


6.6 


-6.340 


0.020 


-0.007 


[20,24] 


1 


^.619 


0.007 


-0.005 


[26] 


1.5 


-3.073 


0.005 


-0.004 


[20,23,25]* 


8 


-2.066 


0.002 


-0.002 


[23]* 


9 


-0.989 


0.001 


-0.002 


[25] 


10 


-0.5130 


0.0002 


-0.001 


[23]* 


12 


-0.166 


0.0010 


0.000 


[25] 


15 


-0.0423 


0.0002 


0.000 


[25] 



evaluated only at r = 5.6 bohr. Without additional 
results at other r we decided, for consistency to omit 
these corrections from the determination of 0o7- The 
sum of these corrections is small compared with the 
scatter of the r = 5.6 bohr potentials in Fig. 3, but of the 
same order as the assigned uncertainty 

The model potential defined by Eqs. (l)-(3) was fit 
to the sum of two quantities, the selected potentials and 
the corresponding DBOC. The input potentials were 
weighted by the inverse squares of the uncertainties 
t/(0) in the fit. The coefficients determined in the fit are 
listed in Table 2. The variance of the fit residuals in the 
determination of 0o7 was 0.6. 

The upper part of Fig. 1 shows the potential 0o7 and 
the selected data used in its determination. The lower 
part of Fig. 1 and Fig. 3 show fractional differences 
between many recent theoretical potentials and 0o7- 
Figure 2 shows the normalized residuals {<p- <p^)IU{<p). 

To assess the uncertainty of 0o7 and the propagation 
of this uncertainty into computed thermophysical prop- 
erties, the potentials 0o7+ and 0o7- were developed. The 
potential was refitted to theoretical potentials shifted by 
their uncertainties, that is, to + A^dbqc ^ U{<p). 
Similarly, the effects of the diagonal Born- 
Oppenheimer correction were assessed by determining 



the potential ^^boc through fits to the theoretical values 
without adding the correction. 

The uncertainty of 0o7 is difficult to quantify Figure 
2 shows that all but one of the theoretical potentials 
used in fitting 0o7 differs from 0o7 by less than the cor- 
responding uncertainty, consistent with the fit variance 
of 0.6. Figure 3 shows that all of the theoretical values 
at r = 4 bohr and r = 5.6 bohr that were used in the fit 
either fall in the range between 0o7- and 0o7+ or have 
uncertainties overlapping this range. These observa- 
tions suggest that the uncertainty in 0o7 should be inter- 
preted as having a large coverage factor [29] ^„^2. 
Table 4 summarizes the properties of the potentials 
used in this work, and the bound state energies (for 
angular momentum index ^ = 0) determined from the 
potential. 



3, Atomic Interactions 

The thermophysical properties of helium can be 
evaluated using the formalism of quantum statistical 
mechanics. In particular, the virial coefficient of the 
equation of state, the viscosity, and the thermal conduc- 
tivity can be expressed in terms of the phase shifts 
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7^ (bohr) 



Fig. 1. Top: The model potential 0o7 (solid line) and theoretical values of (open circles) used in its determination. 
(The vertical scale is proportional to sinh" (50/K), which is approximately logarithmic for large and linear for small 
|0|.) Bottom: Fractional differences between theoretical values of and the model potential 007^ with error bars as 
assigned by the authors (when available). The data sources are D [17], * [18], x [19], + [20], A [23], O [24], T 
[25], V [26]. The potentials 0o7± and 0^)0 [1] are shown as solid lines. 
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Fig. 2. Differences between the theoretical potential values used in fitting 0o7 and the potential, divided 
by the uncertainties of the potential values (See Table 3). 
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Fig. 3. Comparisons of theoretical values of with the (unretarded) model potential (^7. The 
plotted values, arranged chronologically by date of publication, are from the following sources: 
* [18], X [19] (upper bound), + [20], D [21], A [22], A [23], T [25], V [26]. The dotted lines 
represent 0q7±. These bounds encompass the eight values of 0(r) published since 1999 
[19-23,25,26], or overlap the authors' ^„= 1 uncertainty estimates. The values of 0o7 less the 
diagonal Bom-Oppenheimer correction are (292.64 ± 0.13) K at 4 bohr and (-10.996 ± 0.015) K 
at 5.6 bohr. 
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Table 4. The potential minima (1)^^^ = (pir^) for the potentials used in this work, and the corresponding 
bound-state energies. (The retardation corrections^ were included in these calculations.) 



Potential 


0min(K) 


fm (bohr) 


He mass 


Abound (mK) 


<^07 


-10.999 


5.608 


Atomic 


-1.555 


<^07- 


-10.983 


5.608 


Atomic 


-1.667 


<^07+ 


-11.014 


5.607 


Atomic 


-1.438 


<^nboc 


-10.985 


5.608 


Atomic 


-1.550 


^01 


-10.999 


5.608 


Nuclear 


-1.520 



associated with the interaction of a pair of hehum 
atoms. The theory and equations used in determining 
the thermophysical properties are summarized in Sec. 
3.1. The following section 3.2 describes the computa- 
tional techniques used to determine the thermophysical 
properties. 



X^ = kt4 [cos d^ ■ j] (Kr) -sin d^ ■ y^i Kr)l (7) 

where7^(^ and>'^(^ are spherical Bessel and Neumann 
functions. For large Kr the asymptotic form of ;t^^. is 



Xi ,^^> A^sm(Kr-in/2-^d^X 



(8) 



3.1 Formalism 

The interaction of two atoms with a spherically sym- 
metric potential 0(r) is described by a quantum 
mechanical wave function ^'^(r)Y^^/r, where r is the 
separation distance and Y^^ is a spherical harmonic. The 
radial function satisfies 



dr' 



-^^-^[m-E]\^,(r)=0, (4) 



where jU is the reduced mass of the He-He system, m^ is 
the electron mass, lengths are expressed in units of the 
Bohr radius ag, and energies are expressed in units of 
hartree (Ej^). 

The solutions to Eq. (4) fall into three ranges. For 
small r, where the potential is much larger than the 
angular-momentum term, the solutions must be deter- 
mined numerically In the second region of intermedi- 
ate r, the potential is negligible but the angular momen- 
tum term is significant, so Eq. (4) takes the form 



dr' 



^(^+1) |,,2 



X,ir)=0, 



(5) 



where 



K'=i2n/m^)E, 



(6) 



that is, K is just the wave number k = kIqq in atomic 
units. The general solution of Eq. (5) is 



which can be recognized as the solution to Eq. (4) in the 
third region, where both the potential and the angular 
momentum term are negligible. The thermophysical 
properties of interest depend on the phase shifts 5i{E), 
The virial coefficient of "^He depends on the sum 



S{K)= ^ {21 + \)5,{K). 

^=0,2,4,... 



(9) 



The convergence of this sum is discussed in the next 
section. The viscosity and thermal conductivity depend 
on the quantum cross-sections [30] which are expressed 
in terms of much more rapidly converging sums. 

3.2 Numerical Techniques 

Numerical solutions of the radial equation (4) were 
determined with Numerov's method using an integra- 
tion step size 



h^=2x\0-' ■E-'^\ 



(10) 



This step size was determined empirically to insure that 
phase shifts obtained with step sizes h^2 or h^A did not 
differ from those determined with step size h^ within 
the error criterion |A5f| < 10"^. Calculations were made 
for a series of discrete energies in the range 10"^^ < 
EIE^ < 1 . The discrete energies were distributed uni- 
formly on a logarithmic scale. 

For each discrete energy, Eq. (4) was integrated 
upward in r, for ^ = 0, 2, . . . iy. A series of nodes of 
^^ (r) were found at coordinates r„, n = 1, 2, . . . . The 
phase shifts at node n, d^^^, defined by 
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i^hn=Jl(^''n)/ytiK''n)- 



(11) 



were determined successively. The asymptotic phase 
shift as r„ — > oo was obtained when the phase shifts 
evaluated at a series of nodes agree to within the preset 
convergence criterion. Convergence was accelerated by 
using the semi-classical (JWKB) approximation 
[31,32]. The convergence criterion was that the stan- 
dard deviation of three successive values of 5^ „ was 
less than 1 0~^. The maximum angular momentum index 
iy was the minimum of either 1000 or the index when 
\di\ became less than 10"^. 

Equation (11) only determines the phase shift within 
an additive multiple of n. Two conditions were used to 
get the total phase shifts needed in the sum (9). (1) The 
limiting values were lim£_^5o(£) = n and lim£_^Q5^ (E) = 
for ^ > 0; and (2) S^ (E) is a continuous ftinction ofE 
[33]. 

Figure 4 shows the dependence of the phase shifts on 
£ and E. It is clear that for small E, the sum (9) is dom- 
inated by the £ = term. For larger E many terms con- 
tribute to the sum. The Bom approximation 



^Bi =■ 



IjlK 



Jo>W[^"^( 



Kr)fr'dr 



(12) 



(see, e.g. Eq. (38.14) of Ref [34]) for the phase shift is 
useful in considering the rate of convergence. For small 
Kr the spherical Bessel ftinction can be approximated 
by the leading term in the Taylor series, {KrfUli + 1)! !. 
The contribution to the integral in Eq. (12) for small r 
thus decreases rapidly with I, The spherical Bessel 
ftinction has a maximum for Kr near 1+ \. For larger l 
the Bom approximation thus becomes dependent main- 
ly on the weaker attractive part of <^r) The contribu- 
tions from power-law potential terms have a simple 
form: 



h.^\liJA^r)fr'^dr, 



which has the values 



3kk' 



and 



{21 - 3)i2i -\)(2i +\)(2i +3)(2-£ +5) 



(13) 



(14) 



^a ~ 



Ak" 



\5{i-2){i-\)i{i+\){i +2){i +3) 



(15) 



100 - 



— 1°"^ 

- 10-4 



10- 



"10 



-8 



1 - 



T 1 1 r 



T r 



T r 









.11 J-i ■■•' I '' \L VL ii 




IQ-IO IQ-S lQ-6 ^Q-4 

Lj (Ha) 

Fig. 4. Representative phase shifts as functions of energy. The phase shifts are all positive for small E\ 
^ has a zero-^ hmit of 7C, otherwise 5/0) = 0. The lines represent, from left to right, ^ = 0, 2, 4, 10, 20, 
40, 100, 200, 400, and 1000. 
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for the dipole-dipole interaction with and without retar- 
dation. The infinite sums of these terms are 



!/« = 



I/n=- 



2,nK' 



4fj(f,+l)(f,+2)(A+3) 



I even 



(16) 



(17) 



These can be used to get upper and lower limits for the 
contributions of the C^fjr^ term to the truncation error 
of the sum (9). The following test was made to check 
the summation error. For £"> 0.001 hartree, Eq. (12) 
was used to obtain phase shifts for t^ = 1000 < £ < ^2 = 
3000, and the corresponding contributions to the sum 
(9) were evaluated numerically. Equations (16) and 



(17) were then used with 



->- 



, to estimate the trunca- 



tion error of these numerical sums. The results so 
obtained were then compared directly with upper and 
lower limits based on Eqs. (16) and (17). The numeri- 
cal sums were found to lie very close to the product of 
Q and the right-hand side of Eq. (16). The reason is 
that asymptotic phase shifts for i = 1000 are obtained 
when r is some multiple of 27i;/k* beyond the first zero 
of the spherical Bessel function y\ooo('^^)? which occurs 
nearr= 1000.5/ k*- 11.7/V£'. For £"> 0.001 hartree this 
is reached before retardation is significant. For smaller 
E, the nodes r„ where Eq. (11) is evaluated occur at 
larger radii where retardation may be important, but the 
phase shifts decline sufficiently rapidly with increasing 
i that convergence is obtained for i <s: 1000. 



4, Virial Coefficients 

The second virial coefficient of "^He is [33] 

^ = ^th+^ideal+^ound' 



where 



B^=-2N,K'aIJ{KT), 



B,,^=-NA'I\6, 



B^.,=-NA'W^"-^l 



Jo 



(18) 

(19) 
(20) 
(21) 
(22) 



0^ = imJm^)iEjKj,). 

In these equations, A = V2Ar , where 

h 



K 



4'^m^kJ 



(23) 



(24) 



is the thermal de Broglie wavelength, and -JJ, is the 
bound state energy in K (Table 4). The temperature 
derivatives of 5(7) can be evaluated directly from Eqs. 
(18)-(24). Numerical evaluation of the derivatives 
requires the integrals I2 and I^ in addition to /q. 

The thermal contributions B^{T) require numerical 
integration over k. The integrals could formally be 
written with E as the integration variable. However, the 
dependence of the sum terms in the integrand for small 
K was found to be approximately linear in k* oc ^^ so a 
better spline approximation was obtained by using k as 
the independent variable. 

Formally, the upper limit of integration is infinite. In 
practice, the phase shifts become increasingly difficult 
to calculate at higher energies. Calculations were made 
only up to £" = 1 hartree. The argument of the exponen- 
tial factor in the integrand, -a/c^/T, has a maximum 
value at £" = 1 hartree equal to -3.16 x 10^ YJT, so the 
integrand is vanishingly small at k^^, even at T = 
10000 K (exp(-31.6) - 1.9 x 10"^'). The upper limit of 
integration can thus be safely set at k,.^. 

Numerical integrations were required for the inte- 
grals /q, h, and 74. For each case, the sum S(k*) was 
approximated by cubic splines. The number of knots 
per decade of energy E was 40 for all except £" > 0. 1, 
where 80 knots were required in order to resolve the 
rapid dependence of the phase shifts. The integrals 
were calculated as the sum of a series of integrals with 
?c-limits 0-0.01, 0.01-0.1, 0.1-1, 1-10, and lO-K*^. 
This procedure insures sufficient sampling of the inte- 
grands, whose peak values depend strongly on T. 

Figures 5-7 and Table 5 show the virials and the first 
two temperature derivatives as calculated in this work. 
Note that the effects of <^07± on the results is approxi- 
mately symmetrical. Half of the difference of each cal- 
culated property, as computed with <^07+ and <^07-? was 
chosen as a conservative {k^ ~ 2) estimate of the uncer- 
tainty U{x) of property x. These uncertainty estimates 
are well-represented by functions of the form 



kJJ{x) =1 cm^mol ^ exp 



Y^cXHTi K)r 



L«=o 



(25) 



and 
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Fig. 5. Virial coefficient of He, calculated under various assumptions. The plots of ±5(7) = ±Bqj(T) were calculated with ^7 and 
atomic masses. The plotted differences are A^ = B^- Bqj, where x designates the way the virials were calculated; x = 00, 07±, and nboc 
indicates the use of atomic masses and the potentials 0qq, 0q7±, and 0nbocJ ^ ^ ^^"^ indicates calculations with 0q7 and nuclear, rather than 
atomic masses. 
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Table 5. Thermophysical properties of He calculated in this work. Calculated quantities are printed with at least one 
excess figure as an aid in smooth interpolation; for the uncertainties ofB and its derivatives use Eqs. (25) and Table 6; 
for the uncertainties of r] and A use Eq. (37). 



T 


B 


TB' 


f^u 


V 


A 


(K) 


(cm moF ) 


(cm^mol"^) 


(cm^moF^) 


(juPa-s) 


(mWm"^K"^) 


1.0 


^75.05 


669.19 


-1790.42 


0.3292 


2.632 


1.2 


-369.75 


495.66 


-1294.92 


0.3405 


2.720 


1.4 


-301.99 


388.71 


-986.65 


0.3583 


2.845 


1.6 


-254.99 


318.11 


-783.24 


0.3844 


3.033 


1.8 


-220.53 


268.92 


-642.77 


0.4183 


3.283 


2.0 


-194.14 


233.07 


-542.07 


0.4586 


3.588 


2.25 


-168.70 


200.19 


^51.94 


0.5161 


4.030 


2.5 


-148.92 


175.86 


-387.37 


0.5791 


4.519 


2.75 


-133.07 


157.15 


-339.36 


0.6457 


5.037 


3.0 


-120.06 


142.28 


-302.48 


0.7141 


5.569 


3.5 


-99.90 


120.06 


-249.71 


0.8511 


6.637 


4.0 


-84.96 


104.14 


-213.73 


0.9834 


7.666 


4.5 


-73.42 


92.10 


-187.46 


1.1078 


8.636 


5.0 


-64.23 


82.63 


-167.31 


1.2239 


9.542 


6.0 


-50.48 


68.63 


-138.19 


1.4339 


11.184 


7.0 


^0.683 


58.72 


-117.98 


1.6209 


12.650 


8.0 


-33.346 


51.328 


-103.05 


1.7919 


13.992 


9.0 


-27.646 


45.582 


-91.55 


1.9514 


15.244 


10 


-23.090 


40.984 


-82.39 


2.1023 


16.427 


11 


-19.366 


37.216 


-74.93 


2.2463 


17.556 


12 


-16.267 


34.069 


-68.73 


2.3846 


18.641 


14 


-11.407 


29.105 


-58.988 


2.6472 


20.699 


16 


-7.776 


25.358 


-51.679 


2.8947 


22.638 


18 


^.965 


22.423 


^5.981 


3.1299 


24.481 


20 


-2.729 


20.060 


^1.408 


3.3550 


26.244 


22 


-0.911 


18.113 


-37.653 


3.5715 


27.939 


23 


-0.125 


17.262 


-36.014 


3.6768 


28.764 


24 


0.592 


16.479 


-34.509 


3.7804 


29.575 


25 


1.250 


15.757 


-33.121 


3.8824 


30.374 


26 


1.855 


15.088 


-31.837 


3.9828 


31.160 


28 


2.928 


13.888 


-29.537 


4.1795 


32.700 


30 


3.850 


12.841 


-27.534 


4.3709 


34.199 


35 


5.663 


10.729 


-23.496 


4.8301 


37.793 


40 


6.986 


9.123 


-20.432 


5.2659 


41.204 


45 


7.985 


7.860 


-18.024 


5.6828 


44.466 


50 


8.758 


6.838 


-16.077 


6.0837 


47.603 


60 


9.860 


5.286 


-13.117 


6.8465 


53.570 


70 


10.586 


4.1591 


-10.967 


7.5673 


59.208 


80 


11.0827 


3.3033 


-9.329 


8.2549 


64.585 


90 


11.4314 


2.6306 


-8.038 


8.9149 


69.746 


100 


11.6795 


2.0876 


-6.994 


9.5519 


74.726 


120 


11.9830 


1.2650 


-5.403 


10.7689 


84.240 


140 


12.1311 


0.6715 


^.2464 


11.9245 


93.272 


160 


12.1903 


0.2235 


-3.3667 


13.0310 


101.919 


180 


12.1956 


-0.1262 


-2.6743 


14.0968 


110.248 


200 


12.1673 


-0.4063 


-2.1150 


15.1284 


118.308 


225 


12.1026 


-0.6863 


-1.5506 


16.3769 


128.062 


250 


12.0183 


-0.9096 


-1.0953 


17.5862 


137.510 


273.16 


11.9301 


-1.0791 


-0.7458 


18.6765 


146.027 


275 


11.9228 


-1.0913 


-0.7205 


18.7620 


146.695 


298.15 


11.8289 


-1.2315 


-0.4280 


19.8245 


154.994 


300 


11.8212 


-1.2418 


-0.4065 


19.9084 


155.649 


325 


11.7167 


-1.3680 


-0.1398 


21.0288 


164.400 


350 


11.6113 


-1.4752 


0.0893 


22.1260 


172.970 


375 


11.5063 


-1.5671 


0.2883 


23.2024 


181.375 


400 


11.4026 


-1.6465 


0.4626 


24.2598 


189.633 
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Table 5. Thermophysical properties of He calculated in this work. Calculated quantities are printed with at least one 
excess figure as an aid in smooth interpolation; for the uncertainties ofB and its derivatives use Eqs. (25) and Table 6; 
for the uncertainties of r] and A use Eq. (37) — continued. 



T 


B 


TB' 


f^u 


V 


A 


(K) 


(cm moF ) 


(cm^mol"^) 


(cm^moF^) 


(juPa-s) 


(mWm"^K"^) 


450 


11.2008 


-1.7763 


0.7531 


26.3241 


205.753 


500 


11.0082 


-1.8771 


0.9850 


28.3298 


221.413 


600 


10.6523 


-2.0209 


1.3306 


32.1959 


251.596 


700 


10.3332 


-2.1157 


1.5740 


35.9049 


280.550 


800 


10.0462 


-2.1802 


1.7529 


39.4880 


308.517 


900 


9.7867 


-2.2251 


1.8887 


42.9671 


335.670 


1000 


9.5505 


-2.25660 


1.9943 


46.3584 


362.135 


1200 


9.1354 


-2.29376 


2.1455 


52.924 


413.367 


1400 


8.7804 


-2.31003 


2.2454 


59.256 


462.771 


1600 


8.4715 


-2.31430 


2.3139 


65.402 


510.71 


1800 


8.1990 


-2.31131 


2.3617 


71.395 


557.46 


2000 


7.9559 


-2.30378 


2.39552 


77.259 


603.20 


2500 


7.4448 


-2.27457 


2.44221 


91.470 


714.02 


3000 


7.03314 


-2.23916 


2.45852 


105.185 


820.95 


3500 


6.69073 


-2.20239 


2.45921 


118.527 


924.97 


4000 


6.39902 


-2.16616 


2.45129 


131.578 


1026.70 


4500 


6.14591 


-2.13125 


2.43847 


144.394 


1126.60 


5000 


5.92310 


-2.09794 


2.42281 


157.018 


1224.99 


6000 


5.54615 


-2.03623 


2.38740 


181.810 


1418.18 


7000 


5.23652 


-1.98063 


2.35026 


206.135 


1607.72 


8000 


4.97538 


-1.93035 


2.31346 


230.116 


1794.54 


9000 


4.75070 


-1.88461 


2.27786 


253.835 


1979.31 


10000 


4.55433 


-1.84276 


2.24378 


277.355 


2162.52 



with coefficients listed in Table 6. The table also 
includes an uncertainty calculation for the acoustic vir- 
ial 

P, = 25 + 2(r„ -l)TB' + iro -\fT'B'l y„ (26) 

where y^ = 5/3 for helium. Equation (25) represents the 
uncertainties ofB, TB\ and T^B'' within 2 %, 3 %, and 
2 % (rms), respectively, and with a maximum error less 
than 10 %. The uncertainty of ^3^ is represented within 
2 % (rms) with a maximum error of 5 %. As noted pre- 
viously, the uncertainty of 0o7 has a large coverage fac- 
tor K^^l', a similar coverage factor apphes to the 
uncertainties expressed in Eq. (25) and Table 6. 

Figures 5-7 show that that the effects of neglecting 
the diagonal Bom-Oppenheimer correction are no larg- 
er than the uncertainties so assigned, and that the effect 
of using nuclear rather than atomic masses is less than 
the uncertainties except at the highest temperatures. 
The differences between values of 5(7) calculated with 
007 and 000 [1] differ by less than the combined uncer- 
tainties (Eq. (25), Table 6 of Ref [1]). 

As noted above, the integrals required for B^ and its 
temperature derivatives were done numerically. The 
automatic integration routine was controlled by speci- 



fying an error criterion, which was set sufficiently low 
that errors in B and its derivatives due to the numerical 
integrations were negligible. This process only insures 
that the spline-approximated integrand is integrated 
accurately It is of more interest to insure that the 
approximation of the sum term by a spline does not 
introduce a significant error into the calculation. To 
estimate this, the number of knots was reduced by elim- 
inating alternate knots, and recalculating B{T) and its 
derivatives with the cruder spline approximation. The 
absolute fractional differences of the two evaluations of 
B{T) was less than 3x10"^ except near the zero of 5(7). 
The maximum absolute fractional differences of the 
two evaluations of TB\T) was 4 x 10"^, and the maxi- 
mum absolute fractional differences of the two evalua- 
tions of r'5''(7) was 2 x 10^. 

We recommend the use of cubic spline interpolation 
for estimation of B{T) between temperatures listed in 
Table 5. Our tests of such interpolations indicate that 
the fractional interpolation error is generally much less 
than 10^ except at the temperature extremes. 
(Interpolation near the extremes can be improved by 
using the tabulated higher derivatives to set the end 
conditions.) 
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Table 6, Coefficients in Eq. (25) for estimating the uncertainty of 5(7) and its temperature derivatives. 
Property Cq Cy c^ c^ c^ 



B 


0.1341 


-1.4474 


0.0960 


-0.00327 


- 


TB' 


0.6612 


-1.8415 


0.2173 


-0.02476 


0.00128 


f-B'' 


1.8238 


-2.2109 


0.3379 


-0.04263 


0.002166 


jS« 


0.2661 


-1.4560 


0.1134 


-0.00479 


- 



5, Viscosity and Thermal Conductivity 

The kinetic coefficients depend on the quantum 
cross-sections [30] defined by 



Q 



,(2) 



%n^{2l + \){i+2) 



1^ i=o 



ti 2^ + 3 

2(i + \)(i + 2)(2f+6i-3) 



sin'(5, -5,^^), (27) 



(2i-\)(2i+3)i2i+7) 

(^ + l)(^ + 2)(^+3)(£+4) 



sin' (5, -5,^2) 



(2i + 3)i2i+5)i2i+7) 



sin (d, -5,^,) 



(28) 



and 






15(^+l)(^ + 2)(r +6r+t -24^+9) 
_(2^-3)(2^-l)(2^+3)(2^+7)(2^+9) 



sin (3,-3^^^) 



+ sm (6,-6,^^) 



{2£-i){2£ + 3){2i+5){2i+7){2i+n) 
(i + l)(l + 2)(l+3)(i+4)(i+5)(i+6) 
i2l + 3)i2i + 5)i2i + 7)i2l+9)i2i+n) 



sin (6^-6^^^) 



. (29) 



Equations (27)-(29) converge rapidly; numerical 
evaluation was straightforward. The collision integrals 
needed for computation of kinetic coefficients are 
expressed in terms of normalized cross sections, 
defined for even n>Oby 



-)(«)* 



-,(n) 



nr^^nl{n + \) 



(30) 



where r„ (actually an arbitrary length) is the radial posi- 
tion of the potential minimum. The collision integrals 
are defined as 



Q 






(31) 



where l3 = E}^/(kBT). Colhsion integrals with n = 2, s = 
2, 4, ... 10; 77 = 4; 5 = 4, 6, 8; and n = 6, ^ = 6 are need- 
ed for the fifth-order calculation of viscosity and ther- 



mal conductivity [35]. The collision integrals were cal- 
culated by using cubic spline representations of the col- 
lision integrals, and dividing the integrals in Eq. (31) 
into 11 sub-intervals, with limits 0-1 0"^^ 10"^^-10"^ 
. . . 0.1-1. This division insured adequate sampling of 
the integrands, whose peak locations vary rapidly with 
temperature. (The errors introduced by truncating the 
infinite integral are neglibible.) 
The viscosity is [35,36] 



7] 






(32) 



where /^^"^ is obtained by solving a set of linear equa- 
tions 



^^11 hi hs 



H = 



^4 



hs\ 


rn 




(\\ 


h 


4 







hs 


4 


= 





h 


^4 







hs 


4 1 




,0, 



= e, 



(33) 



for ^y =f^"^lbii. The components of the symmetric 
matrix B are listed in Appendix A of Ref [35]. In par- 
ticular, since b^i = AOP'^^^, the viscosity can be 
expressed as 



S^nmk^T 



4;rr: 



(34) 



Similarly, the thermal conductivity can be determined 
from the solution of 



AC=e„ 



(35) 



where the components of A are defined in Appendix B 
of Ref [35], and f is a column vector with components 
Q. The thermal conductivity depends only on f^: 



_15k^^nmk^T 
A- — — Qy. 



I6mr^ 



(36) 
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To insure that the complicated formulas for the com- 
ponents of B (and the corresponding matrix for the 
thermal conductivity) were transcribed accurately the 
following procedure was followed. The formulas were 
extracted from an electronic copy of Ref [35]. These 
were further edited to conform with Fortran notation. 
Subsequently Viehland provided Fortran codes that 
generated the Fortran code for calculating the matrices 
directly [37]. Numerical evaluations using the two 
implementations were identical within machine preci- 
sion. 

Errors in the numerical integrations required for cal- 
culating rj and X were estimated by eliminating alter- 
nate knots in the spline representations of the collision 
integrals and repeating the calculations. The two values 
of T}(T), and the two values of A{T) so determined had 
an absolute fractional difference of less than 3 x 1 0"^ at 
r = 1 K. This difference declined with T and remained 
below 1.2 X 10-^forr>20K. 

The viscosities and thermal conductivities deter- 
mined in this work are listed in Table 5. Figure 8 and a 
nearly identical figure for AA /A show the sensitivity of 
the calculations to the choice of potential and the 
choice of nuclear instead of atomic masses. The effects 
of using potential 0o7± is nearly symmetric. Half of the 
differences between values of rj or X calculated with 
these two potentials approaches 0.35 % at low temper- 
ature. The differences reverse sign near 42 K. Above 
this temperature, the half-difference is bound by 
0.02 %. A reasonable estimate of the the relative uncer- 
tainty U^ in either rj or X is the minimum of 0.35 % and 
the equation 

k,UXv) = KUIX) = 0.0002 + 0.005 KIT. (37) 

Values of the viscosity and thermal conductivity at 
temperatures between those listed in Table 5 can be 
obtained by interpolating with cubic splines. Our tests 
indicate that cubic spline interpolation introduces a 
fractional error of less than 10"^ except near the temper- 
ature extremes. 



6, Validation of Computations 

The Fortran code used for calculating the phase 
shifts and for subsequent calculation of the thermo- 
physical properties was tested by an independent devel- 
opment of new codes by one of us (Mehl) to test the 
results of Hurly and Moldover [1]. The test demonstrat- 
ed excellent agreement of the sum (9) and the quantum 
cross-sections (27)-(29). 



The test revealed two errors in the calculation of the 
thermophysical properties reported in Ref [1]. The first 
was an incorrect sign assigned to the bound-state con- 
tribution to the published virials, which mainly affect- 
ed the low temperature results for "^He and for ^He-'^He 
mixtures. The second was due to inconsistent units con- 
version. The code used by Hurly to calculate the ther- 
mal conductivity was based on the equivalent of Eq. 
(36) in Hirschfelder et al. [36]. Their Eq. (8.2-31) uses 
a calorie unit in a numerical prefactor. Conversion of 
this to J using a current definition of the calorie intro- 
duced a factor of 1.000545 error in the thermal conduc- 
tivity results pubhshed in Ref [1]. The published val- 
ues are high by this factor. 



7, Comparisons With Recent 
Experiments 

Hurly and Moldover [ 1 ] compared their results with 
a wide range of experimental results. Here we limit our 
comparisons to a few accurate experiments published 
since 2000. Figure 9 compares the recent second virial 
measurements of McLinden and Losch-Will [38]. The 
agreement is excellent. 

Figure 1 compares the recent measurements of the 
acoustic virial by Pitre, Moldover and Tew [5]. The 
measurements fall well within the combined (k^ = 2) 
uncertainty of the predicted slope P^ and the experi- 
mental uncertainty except at high temperatures, where 
the disagreement is on the order of the scatter in the 
measurements. 

Berg's high quality measurement of the viscosity 
[39,40] at 298. 15 K (expressed with a ^„ = 2 uncertain- 
ty), (19.842 ±0.014) jUPa-s, and the calculated value 
(19.824 ± 0.004) jUPa-s differ by the sum of their k, = 2 
uncertainties. 



8, Concluding Remarks 

As shown in Fig. 3, multiple research groups have 
provided us with very accurate ab initio "data" at 4.0 
and 5.6 bohr. In order to fully exploit these data, it 
would be desirable to have theoretical potentials of 
similar accuracy at nearby r. The most demanding gas 
metrology is conducted near 273 K; thus, it would be 
very desirable to generate ab initio values of the poten- 
tial at, for example r = 3.89 and 4.13 bohr (correspon- 
ding to = 200 K and 450 K) with uncertainties com- 
parable to those already achieved at 4.0 and 5.6 bohr. 
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Fig. 8. Sensitivity of tlie viscosity of He to various options in the calculations. The fractional difference between ?7^ and the value cal- 
culated with 0Q7 and atomic masses is plotted as the fraction Arf/ri = {t}^- riQj)(ri^, where x specifies the type of calculation (See cap- 
tion to Fig. 5.) A similar plot for the thermal conductivity differs from this plot only in minor details. 
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indicate the experimental ^„ = 1 uncertainties. 
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Fig. 10. Acoustic virial coefficient of Pitre et al. [5] compared with values calculated with 0q7; Ap^ = 
Pa,cxpt~ PaMc ^^^ dashed lines are plots of P^qj^- P^qj, and indicate the uncertainty of the theoretical 
calculation. The error bars indicate the experimental (k^ = 1) uncertainties. Other lines show Ap^ corre- 
sponding to 0nn^, and 0ni^. The acoustic virial is clearly sensitive to the differences between the various 
potentials. 
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